Mott insulators of ultracold fermionic alkaline earth atoms in three dimensions 
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We study a class of SU{N) Heisenberg models, describing Mott insulators of fermionic ultra-cold 
alkaline earth atoms on the three-dimensional simple cubic lattice. Based on an earlier semiclassical 
analysis, magnetic order is unlikely, and we focus instead on a solvable large-A'^ limit designed to 
address the competition among non-magnetic ground states. We find a rich phase diagram as a 
function of the filling parameter k, composed of a variety of ground states spontaneously breaking 
lattice symmetries, and in some cases also time reversal symmetry. One particularly striking example 
is a state spontaneously breaking lattice rotation symmetry, where the cubic lattice breaks up into 
bilayers, each of which forms a two-dimensional chiral spin liquid state. 



I. INTRODUCTION 

Ultracold atom experiment techniques enable us to 
vary parameters of quantum many-body systems that 
can hardly be changed in solid state materials.^"' For 
example, in solid state systems the crystal structure is 
selected by nature, so it is usually not easy to study 
the dependence of the system properties on the lattice 
structure. But in ultracold atom experiments the opti- 
cal lattice can be chosen artificially, and its dimension 
and geometry can be varied. Also, we have significant 
freedom to select the constituent particles of a many- 
body system. They can be atoms or molecules, bosons 
or fermions, and so on. Different atoms or molecules in- 
teract with one another quite differently, and in some 
cases the interactions can be tuned with electric or mag- 
netic field. So cold atoms promise to allow us explore 
systems in new parameter regimes, or even systems that 
have no analog in solid state materials. 

Fermionic^ ultracold alkaline earth atoms (AEAs) have 
attracted significant interest recently due to their unique 
properties,''"-^ ' and experimental progress developing the 
study of many-body physics in AEA systems has been 
rapid^''"^^. One key feature of AEAs is the presence, to 
an excellent approximation, of SU{N) spin rotation sym- 
metry, where N = 21+1 and / is the nuclear spin."''"' This 
occurs in both the ^Sq ground state and a metastable ^Pq 
excited state, where the electronic angular momentum 
Je = and the hyperfine interaction is thus quenched. 
This leads to the nuclear-spin-independence of the s- 
wave scattering lengths between AEAs, and to SU{N) 
spin rotation symmetry. When loaded in optical lattices, 
AEA systems are described by S'f/(A^)-symmetric Hub- 
bard models."' Since the largest / obtained using AEA is 
J = 9/2 in the case of ^'^Sr, TV < 10 is the experimentally 
accessible regime. Different setups are possible, and as 
a result, SU{N) versions of several models, such as the 
Kugcl-Khomskii model, the Kondo lattice model, and the 
Heisenberg spin model, can be realized with AEAs as spe- 
cial or limiting situations of the more general Hubbard 
model. 

Among these models, we focus in this paper on SU {N) 
antiferromagnetic Heisenberg models, which describe the 
Mott insulator phase of fermionic AEAs in optical lat- 



tices. More specifically, we are concerned with such mod- 
els on three dimensional lattices, which have received 
much less attention than the one- and two-dimensional 
cases. Because of the enlarged symmetry, the number of 
spins needed to make a singlet, denoted by fc, is in general 
larger than two. In the simplest AEA Heisenberg model 
with one atom per lattice site, k = N. In addition, in 
the semiclassical limit of the Heisenberg models that can 
be realized using AEAs, two neighboring classical spins 
prefer energetically to be orthogonal rather than anti- 
parallel.' Both these features contrast with SU{2) an- 
tiferromagnetic Heisenberg models appropriate for some 
solid state materials, where neighboring pairs of spins 
can and tend to form singlet valence bonds, and neigh- 
boring classical spins prefer to be anti-parallel. We can 
thus expect new physics in SU{N) Heisenberg models 
with fc > 2. 

Indeed, Ref. 7 argued that the underconstrained nature 
of the semiclassical limit makes magnetic order unlikely 
for large enough N on any lattice, and non-magnetic 
ground states are more likely. While the models of phys- 
ical interest are challenging to study directly, informa- 
tion about possible non-magnetic ground states can be 
obtained in a large- A" limit designed to address the com- 
petition among such states."'^""'"' Such a large-A^ study 
was carried out for AEA SU{N) Heisenberg models on 
the two-dimensional square lattice in Refs. 7 and 1.3. One 
possible non-magnetic state is a cluster state, where clus- 
ters of k (or a multiple of fc) neighboring spins form sin- 
glets; this is a generalization of a valence bond state. An- 
other possibility is a spin liquid state, where full trans- 
lational symmetry is preserved. For the simplest AEA 
Mott insulators (with ^Sq ground state atoms only), on 
the square lattice the large- TV study finds cluster states 
for fc < 4, and a chiral spin liquid (CSL) state for 
fc > 5.''^' The CSL spontaneously breaks time-reversal 
(T) and parity (V) symmetries, and can be viewed as 
a magnetic analog of the fractional quantum Hall effect 
(FQHE), with similar exciting properties of quasiparti- 
cles with anyonic statistics, gapless chiral edge states, 
and so on.^'""'' CSLs have also been found in a variety 
of other exactly solvable models. ""^"^"^ 

The CSL is, however, intrinsically a two-dimensional 
phenomenon, so it is natural to ask about non-magnetic 
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ground states of SU{N) antifcrromagnetic Heisenberg 
models in three dimensions. In this paper, we address 
this question by a large- study of a class of SU{N) 
Heisenberg models on the simple cubic lattice, and find 
a rich phase diagram as a function of k including clus- 
ter states, but also more intricate inhomogenous states. 
Most strikingly, for fc = 7, 10 we find a bilayer CSL state, 
where the lattice spontaneously breaks into weakly cou- 
pled square bilayers (thus breaking rotational symmetry), 
each of which is a two-dimensional CSL. We thus find 
that the CSL survives to three dimensions, relying on 
spontaneous symmetry breaking that results in effective 
quasi- two- dimensionality. 

We now define our model before briefly surveying some 
related prior work. We consider a fermionic AEA with N 
spin species, and put m ^Sq ground state atoms on each 
site of a simple cubic lattice (see Sec. II for more details). 
The atoms form a Mott insulator due to repulsive on- 
site interactions. For simplicity, we consider the case of 
dominant on-site interaction, so that the spin degrees are 
governed by a antifcrromagnetic superexchange interac- 
tion restricted to nearest neighbors. While m = 1 is the 
most interesting situation since it best avoids three-body 
losses, we also consider more generally the case where 
^ is an integer. Then, the minimum number of spins 
needed to make a SU{N) sing let is fc = ^. We some- 
times refer to fc as the filling parameter. When m = 1, 
each spin transforms in the fundamental representation 
of SU (N). In the large- iV limit, N is taken large while fc 
is held fixed. Given the physical interpretation of fc, we 
thus view the large- results for a given fc as a guide to 
the physics of the physically realizable model with m ~ 1 
and N ^ k. 

Our focus is on three spatial dimensions, but we note 
that one-dimensional SU (N) Heisenberg spin chains have 
been solved exactly for the case m = 1,'"" and the effec- 
tive field theory of such chains is understood for general 
m."''^ The latter analysis shows that gapless states with 
quasi-long-range order, as well as gapless cluster states, 
occur in one dimension. In two dimensions, early stud- 
ies of SU{N) antiferromagnets focused on models where 
two neighboring spins can be combined to form a singlet. 
This work included the models we consider for the case 
m = A^/2,'-'' ' but also other SU{N) antiferromagnets 
with spins transforming in two distinct conjugate repre- 
sentations on the two sublattices of a bipartite lattice. ' ' 
Models with fc = 2 have also received attention more 
recently, ^"'■^■^■''^''''^ and two dimensional models with fc > 2 
have been studied' (see Ref. l-! for a 
more detailed discussion of some of these prior works). 
The m = 1 . iV = 3 model on the square lattice is magnet- 
ically ordered,^' and there is also evidence for magnetic 
order for m = 1, iV = 4.^- Only a little attention has 
been devoted to the case of three dimensions, but 
we note the high temperature series study of Ref. 65, 
where the m = 1 model on the simple cubic lattice was 
studied for various values of N, and it was found that in- 
creasing N led to a decreased tendency toward magnetic 



order. References 66 and 67 studied effective models for 
four-site singlet clusters on the cubic lattice. Finally, we 
note that high-spin quantum magnets can also be real- 
ized using ultra-cold alkali atoms. While iV-component 
such systems do not generically obey SU (N) spin sym- 
metry, the symmetry is enhanced above S'?7(2),''^ and 
such systems have received significant attention. 

In Sec. II, we review the large- TV solution to our model. 
This is followed by presentation of the large- results for 
fc = 2, . . . , 10 in Sec. Ill, together with a discussion of how 
those results are obtained and checked. As part of that 
discussion, we develop an interesting relation between 
some cubic lattice saddle points (including the ground 
state saddle points for fc = 5, . . . , 10) and saddle points 
on the single-layer square lattice with filling parameter 
fc' = fc/2. The paper concludes with a discussion of the 
striking properties of the bilayer CSL state (Sec. IV). 

II. THEORETICAL MODEL 

The SU{N) Hubbard model 

TiHubbard = 51 i'^r^^r'a + h.C.) 
{rr') 

+ iU/2)Y,icr^Cr^-r>i)\ (1) 

r 

describes the behavior of fermionic AEAs on an optical 
lattice.'' Here c"^ and Cra are the creation and annihila- 
tion operators for the fermionic atom with spin state a at 
site r. The sum in the first term is over nearest- neighbor 
pairs of lattice sites. We will primarily consider the sim- 
ple cubic lattice. We choose the number of atoms so that 
m is the integer number of atoms per lattice site. There 
are N spin states, a,/3 = 1,2,...,A'^, and spin indices 
are summed over when repeated. The total number of 
lattice sites is Ng. The operator c"^ transforms in the 
fundamental representation of SU{N)^ while Cra trans- 
forms in the anti-fundamental representation, which is 
related to the fundamental by complex conjugation. The 
upper and lower positions of the Greek indices are used 
to indicate the distinction between these two representa- 
tions (they are unitarily equivalent only for N = 2). 

As is well known, the SU{2) Heisenberg model can 
be obtained as a low energy effective description of the 
SU{2) Hubbard model when U ^ t. The generaliza- 
tion to the SU{N) version is straightforward. In sec- 
ond order degenerate perturbation theory, one obtains 
the SU{N) antifcrromagnetic Heisenberg model defined 
by the Hamiltonian 

{rr') 

with the Hilbcrt space restricted by f^^fra — fn, and J = 
2t^ /U > 0. We now use /"''' rather than c"''' to denote 
the fermion creation operator, to emphasize that once we 
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pass to the Hcisenberg model, the fermions do not move 
from site to site. This is important, because the structure 
of the large- mean-field theory is that of a hopping 
Hamiltonian for the /"^ fermions, but it is not correct 
to interpret this hopping as motion of atoms. Instead, 
in the large- mean-field theory, the /"^ fermions are 
spinous, fractional particles that may be either confined 
or dcconfincd depending on the nature of fluctuations 
about a mean-field saddle point. See Ref. l.> for further 
discussion of this point. 

On each site, there are m atoms that form a SU{N) 
spin. The Hamiltonian (2) defines an antiferromagnetic 
interaction, since by rearranging the fcrmion operators it 
can be written as 

n^jJ2S^{r)S^{r'), (3) 

(rr') 

where S^{r) = fr^fra flips the spin on site r. 

We study this model on the simple cubic lattice, the 
simplest three dimensional case, with varying parameters 
N and m. While we consider more general parameter 
values, m = 1 is the case of greatest physical interest 
because putting only one atom on each site best avoids 
potential issues due to three body loss. The largest N 
that can be obtained using alkaline earth atoms is iV = 10 
in the case of ^^Sr. 

Based on a scmiclassical analysis, Ref. 7 argued that 
for large enough A^, magnetic ordering is unlikely on any 
lattice. The argument proceeds in the semiclassical limit, 
where a lower bound on the dimension of the ground state 
manifold is derived. For N > Nc, where Nc depends on 
the lattice coordination number, the ground state mani- 
fold is extensive, meaning its dimension is proportional to 
the number of lattice sites. This situation occurs in some 
geometrically frustrated systems and is likely to lead to a 
strong or complete suppression of magnetic order' ', even 
in the semiclassical limit that favors magnetic order by 
construction. Therefore, non-magnetic ground states are 
likely when N > Nc- For the square lattice Nc = 3,' and 
the argument is easily extended to find TVc = 4 on the 
cubic lattice. 

Ideally, we would like to predict the properties of the 
SU{N) antiferromagnetic Heisenberg model on cubic lat- 
tice for < 10, m = 1. But this is extremely challeng- 
ing. Instead, following the work of Refs. 7 and 13, we ap- 
ply a large- A'^ limit in which the model becomes exactly 
solvable, and which allows us to address the competition 
among different non-magetic ground states. We fix the 
ratio fc = ^ (for integer k), while taking both N ^ oo 
and m — >■ oo. We shall sometimes refer to k as the filling 
parameter. For each k we thus obtain a sequence of mod- 
els (A'' = k,m = 1); (A^ = 2k, m = 2), and so on. For 
every model in this sequence, k is the minimum number 
of spins needed to form a singlet, and it is thus reason- 
able that the large- A^ limit may capture the physics of 
the case A^ = fc, m = 1 of greatest interest. 

To proceed with the large- A^ solution, one goes to a 
functional integral representation, where the partition 



function is 

Z = j VxVx'VXDjVfe-^, (4) 

where 

+ / X! [Xrr'Irfr'a+h.C.) 

(r,r ) 

/^A.(/;^/™-m). (5) 

The field Xrr' is a complex Hubbard-Stratonovich field 
that has been used to decouple the exchange interaction, 
and A,, is a real Lagrange-multiplier field enforcing the 
fr^fra = ™ Constraint. The fermion fields / and / are 
the usual Grassmann variables. We have introduced J' = 
NJ; J is held fixed in the large-A^ limit. Finally, = 

Jp^ dr. We shall always be interested in zero temperature, 
i.e. ji — > oo. 

When both A^ and m are large, the effective action 
for X and ^ (obtained upon integrating out fermions), 
is proportional to A^ (since m ~ A"), and therefore the 
saddle point approximation becomes exact for the x and 
A integrals. We can therefore replace x and A by their 
saddle-point values, Xrr' Xrr' and A^ i/ir- The 
saddle-point equations are 

m - {!f!rc) , (6) 

= (7) 

The above averages are taken in the ground state of the 
saddle-point (or mean-field) Hamiltonian 

Hmft=A^E%^+™E^'- 

{rr') r 

+ X! {Xrr'fr^'fr'a + h.C.) - ^ Hrrir, (8) 

{rr') r 

where h,. = /"Vra- 

The ground state is determined by finding the global 
minimum of Emft ({Xrr'} , {/^r}), the ground state en- 
ergy of Haift, as a function of the x's and /u's, with the 
constraint that the saddle point equations must be sat- 
isfied. While any solution of the saddle point equations 
gives an extremum of the energy, in general it is not triv- 
ial to find the global minimum. To address this question, 
we follow Refs. 7 and 13 and apply the combination of 
analytical and numerical techniques developed there, as 
described below in Sec. III. 



III. LARGE-iV GROUND STATES 

In the limit A^ — > oo, the ground states are character- 
ized entirely by the mean-field saddle point values of Xrr' 
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and /ir- The most important information is contained in 
Xrr'j since typically it is possible for a given Xrr' to find 
Hr SO that the density constraint Eq. (6) is satisfied. For 
instance, depending on whether two sites are connected 
(i.e. whether there is a set of nonzero Xrr' 's forming a 
path connecting the two sites), we can tell whether the 
spins on the two sites are correlated or not. Not all the 
information contained in Xrr' is physical. The theory has 
a C/(l) gauge redundancy 

f — f p*<^('") 
Jra r Jra^ 

y , y ,e<'^(O-0('-')) ' 

Xrr' r Xrr'f^ 

so the physical information is contained in the following 
gauge-invariant quantities: (1) magnitude \xrr'\ and (2) 
fiux $ = ai2 + 023 + ^34 + 141 through each plaquette, 
where 1, 2, 3, 4 indicates the four vertices of a plaquette 
and Urr' is the phase of the Xrr', *-e- Xr'r = e^°'^^' |Xrr '|- 

(Since Xr'r = Xrr'^ "-r'r = -ttrr'-) 

Based on a combination of analytical and numerical 
techniques described below, we found the ground state 
configuration of Xrr' and fir for k — 2, . . . , 10. These re- 
sults, which are rigorous for k = 2,3,4, are summarized 
in Table I. Different types of ground states are found de- 
pending on k. In an n-site cluster pattern of Xrr'i the 
lattice is partioned into n-site clusters such that Xrr' 7^ 
only if r, r' lie in the same cluster. We call the corre- 
sponding ground state a n-site cluster state, which can 
be viewed as a generalization of a valence bond state 
(2-site cluster state, in our terminology). Similarly, a bi- 
layer pattern partitions the lattice into bilayers, and Xrr' 
is only nonzero for r, r' in the same bilayer. The corre- 
sponding ground states are called bilayer states. In all 
cases, each bilayer is comprised of two adjacent {100} 
lattice planes. A CSL bilayer is a special kind of bilayer 
state, where in each bilayer 

{X, {rr') lies within either layer; 
^, {rr') connects the two layers. 

Moreover, there is a uniform flux 



through each plaquette lying within the two layers, and 
zero flux through each plaquette perpendicular to the two 
layers. This situation corresponds to a uniform orbital 
magnetic fleld applied perpendicular to the layers. At 
the mean-fleld level, a single CSL bilayer exhibits integer 
quantum Hall effect with i> = 1 for each spin species of 
fra fcrmion. 

To fully understand the different ground states, one has 
to go beyond the iV = oo or mean-field description. At 
the mean-field level, the number of ground state arrange- 
ments of clusters or bilayers on the cubic lattice diverges 
with the system size. For example, there are usually 
many ways to tile the lattice with a given type of n-site 
cluster. Also, in the CSL bilayer state, the direction of 
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2/4-site cluster 


Fig. la 


-0.125 


3 


6-site cluster 


Fig. lb 


-0.0833333 


4 


4/8-site cluster 


Fig. Ic 


-0.0625 


5 


20-site cluster 


Fig. 2a, 2b 


-0.0445021 


6 


12-site cluster 


Fig. 2c, 2d 


-0.0347222 


7 


CSL bilayer 


Fig. 2c, 2f 


-0.0273888 


8 


8-site cluster 


Fig. 2g, 2h 


-0.0234375 


9 


Inhomogeneous bilayer 


Fig. 2i, 2j 


-0.0188265 


10 


CSL bilayer 


Fig. 2e, 2f 


-0.01577 



TABLE I. Ground state saddle-point patterns of Xrr' , and the 
corresponding energies in units of NJNa for = 2, 3, . . . , 10. 
The different types of large- A'^ ground states are described in 
the text, and depicted in figures as indicated. 



fiux can be chosen independently in each bilayer without 
affecting the N = oo ground state energy. Such degen- 
eracies can be resolved by computing the first correction 
(pcrturbative in 1/A'^) to the ground state energy;"'"' these 
calculations are left for future work. 

In cluster states, another important effect of fluctua- 
tions is to confine the fra fermions; the cluster states are 
thus "ordinary" broken symmetry states, without exotic 
excitations. A more extensive discussion of fluctuations 
appears in Ref. 1.'], and the resulting physical properties 
of the CSL bilayer are discussed in Sec. IV. We have not 
considered the effect of fluctuations in the fc = 9 inhomo- 
geneous bilayer ground state. 

We now discuss the mean-field ground states for each 
value of k. We note that, for fc > 5, we cannot rule out 
the possibility that the true ground state is lower in en- 
ergy than the ground state we found. The ground-state 
clusters for fc = 2, 3, 4 are depicted in Fig. 1. These are es- 
sentially the same as found in the two-dimensional square 
lattice,''' ' but going to the three-dimensional cubic lat- 
tice permits a greater variety of clusters for fc = 3, 4. 

It was noted in Ref. 44 that for fc = 2 there is actually a 
continuous family of N — oo ground states, which can be 
seen for a single square plaquette as shown in Fig. la and 
discussed in the figure caption. This continuous ground 
state degeneracy is also resolved by the order- 1/A^ cor- 
rections to the ground state energy."'"' We found that a 
similar continuous degeneracy occurs for fc = 4 on a sin- 
gle cube (see Fig. Ic). As in the figure, consider a single 
cube with flux $t through the top and bottom plaquettes 
{i.e., those lying in the xy-plane), and fiux $s through 
the side plaquettes (i.e., those lying in the xz- and yz- 
planes). Flux passing from the center of the cube to the 
outside is taken positive. In order to reach the ground 
state we must have 2<f>t 4- 4(f>s = ±27r; we choose the 
positive sign without loss of generality. We let = 4u 
and ~ 7r/2 — 2u; a ground is obtained if we restrict 
< w < 7r/2. In this situation the magnitude \xrr'\ will 
generally differ on vertical bonds and other bonds [shaded 
light (pink) and dark (blue), respectively, in Fig. 1]. The 
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energy is minimized and saturates the lower bound when 
= 2Vcosusin7.. (12) 

The ground-state patterns of Xrr' for 5 < fc < 10 arc 
shown in Fig. 2. For k = 5,6,8 we again find clus- 
ter ground states. The case fc = 8 is particularly sim- 
ple; there, each cluster is a fully symmetric cube with 
I Xrr' I constant on every bond, and no flux through the 
cube faces. The fc = 5 and fc = 6 clusters are conve- 
niently thought of as obtained by stacking two single- 
layer clusters vertically, and connecting them via the ver- 
tical bonds. For fc = 5 each cluster is a stack of two 
ten-site T-shaped objects. The fc = 6 clusters are ob- 
tained by stacking two fc = 3 ground state clusters (see 
Fig. lb). In the fc = 5,6 cases, our numerical calcula- 
tions find evidence for a continuous family of degenerate 
ground states within each cluster, as for the 4-site fc = 2 
clusters and 8-site fc = 4 clusters (Fig. 1). Unlike in those 
cases, however, we have not been able to find a simple 
parametrization of the degenerate ground states. 

For fc = 7, 9, 10, we find bilayer ground states, with 
the CSL bilayer saddle point described above occurring 
for fc = 7, 10. The fc = 9 ground state is more com- 
plicated, spontaneously breaking translation symmetry 
within each bilayer. Time reversal symmetry is broken 
as well by a complicated pattern of fluxes. It is interest- 
ing to note that all the 5 < fc < 10 ground states have 
a bilayer structure, as the clusters for fc = 5,6,8 can 
be arranged into bilayers (see right column of Fig. 2). 
In addition, the two square lattice layers of each bilayer 
have identical Xrr', there is zero flux on the "vertical" 
plaquettes connecting the two layers, and the vertical 
bonds have magnitude |xrr'| = J'/k.''' As discussed be- 
low, this simple structure allows us to exploit a useful 
relation with the single-layer square lattice at filling pa- 
rameter fc' ~ fc/2. 

We now describe how the large- ground states were 
determined. As on the square lattice,''' ' the results for 
fc = 2,3,4 are rigorous, and are obtained by applying a 
lower bound on Em ft obtained by Rokhsar for fc = 2,' ' 
and generalized to fc > 2 (with a stronger bound holding 
for bipartite lattices) in Refs. 7 and 13. Cluster states for 
fc = 2,3,4 on the square''^' and cubic lattices saturate 
this lower bound. A necessary condition for saturation 
on a bipartite lattice is that the mean-field single-particle 
energy spectrum must be completely flat, with only three 
energies 0, ±e occuring in the spectrum, and with energy 
— e states filled and others empty. We believe that this 
kind of spectrum can only be produced by a cluster state. 
Moreover, for larger clusters (and thus with increasing fc), 
it becomes harder to arrange for a spectrum containing 
only three energies. While we do not have a rigorous 
proof, we believe saturation is impossible for fc > 4 on 
the square and cubic lattices. 

For fc > 5, we resort to a numerical approach to find 
the ground states. We employ the self-consistent min- 
imization (SCM) algorithm developed in Refs. 7 and 



(a) = 2 




(c) fc = 4 

FIG. 1. Ground-state clusters for k = 2,3,4. Shaded bonds 
are those with Xrr' 7^ 0. Bonds with different shading (or 
color in online version) may have different magnitudes |Xrr'|- 

(a) The fc = 2 ground state clusters are dimers and square pla- 
quettes. The square plaquette is pierced by vr-flux, and the 
ratio of IXrr'l on light (pink online) and dark (blue online) 
bonds can be chosen arbitrarily. Setting |Xrr'| = on the 
two light (pink) bonds breaks the plaquette into two dimers. 

(b) The fc = 3 ground state cluster is a 6-site chain pierced 
by TT-flux. On the cubic lattice, such chains can exist either 
as a flat rectangular loop (left), or as the same loop bent by 
90° in the middle (right). In both cases, Xrr' = on the 
dashed bond passing through the middle of the loop, (c) The 
k — 4 ground state clusters are square plaquettes and 8-site 
cubes with Os-flux through the side plaquettes and <I>t-flux 
through top and bottom plaquettes. There is a continuous 
one-parameter family of ground states on an 8-site cube, de- 
scribed in the text. 

13, which proceeds as follows (see Ref. 13 for more de- 
tails): 

1: Start with = and a randomly generated config- 
uration of Xrr' ■ 
2: Adjust fir to Satisfy the saddle-point equation 

(/r"Vra)=m, for ah r. (13) 

Hr is determined by a multidimensional Newton's 
method.''^ Stop if no solution is found. 
3: Generate a new Xrr' using the saddle-point equation 

Xrr'=-f (/"Vr«). (M) 
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(i) fc = 9 
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FIG. 2. Ground-state saddle point configurations of Xrr' for 
fc = 5, 6, . . . , 10. The right column is a tliree-dimensional view 
of each configuration, with larger magnitude \Xrr' \ indicated 
by darker shading. All these saddle points can be viewed 
as bilayer structures, with Xrr' identical on top and bottom 
layers. The left column thus shows |Xrr'| on a single layer, 
with fluxes indicated except for fc = 9, where the fluxes are 
generally non-zero but follow a complicated pattern. Also, 
for fc = 5, 6 the fluxes and |xrr' I can be changed continuously 
within a single cluster without affecting the energy; only the 
simplest configurations are shown. 



4: Go back to step 2 until Xrr' and /ir converge. 
As long as step 2 is successful, the energy Emft is 
guaranteed to decrease with each iteration of the SCM 
algorithm.''^ ' But a random initial configuration of Xrr' 
does not necessarily converge to the ground state, and 
can instead converge to a local minimum of Emft- 
Therefore, in order to find the ground state, we need to 
try as many independent random initial configurations of 
Xrr' as possible. 

To improve the performance of the SCM algorithm, 
we define Xrr' with fir within some fixed unit cell, 
which is then repeated periodically to cover a finite-size 
Lx X Ly X Lz lattice with periodic boundary conditions. 
For simplicity, we always choose the unit cell to be a rect- 
angular prism with edge lengths lx,y,z (see Fig. 3), with 
primitive Bravais lattice vectors parallel to the edges of 
the rectangular prism.'" For each value of fc, we choose 
the minimum linear system size L = m\n{LxT Ly, Lz) to 
be as large as possible given the constraints of our avail- 
able computing resources and the need to try a reason- 
ably large number of different random initial conditions. 
In some cases we also considered larger system sizes, espe- 
cially when we found competing saddle points very close 
in energy. A more careful study of finite-size effects would 
be desirable, but due to the above constraints we leave 
this for future work. Table II displays the range of unit 
cell dimensions studied for each value of A:, as well as the 
number of random initial conditions tried for each cell, 
and the minimum linear system size L. 

As noted above, the ground states for 5 < ^ < 10 
can all be viewed as bilayer states, which means that 
such saddle points can also be obtained by a studying 
the large-A^ Heisenberg model on a single bilayer. We 
have also carried out SCM numerical calculations in this 
geometry (see Table II and Fig. 3 for more information); 
this is computationally cheaper than the cubic lattice 
SCM calculations, and provides a useful check on those 
results. These bilayer SCM calculations find the same 
ground states as the corresponding cubic lattice calcula- 
tions, except for A: = 9, where the bilayer calculation finds 
a lower-energy state that can then be extended to a cubic 
lattice saddle point. Presumably, this saddle point would 
also be found by SCM on the cubic lattice with enough 
runs using independent random initial conditions. 

There is an interesting relation between certain saddle 
points of a single bilayer, and corresponding saddle points 
of a single-layer square lattice, but with filling parameter 
k' = k/2. The cubic lattice ground states for 5 < fc < 10 
are all of this type. We label the sites of a single bilayer 
by (r, z), where i = 1,2 is the layer index, and r labels 
the square lattice sites within each layer. There are Ng = 
2N'^'^ lattice sites, where TV^'* is the number of sites in a 
single layer. Consider a saddle point where 

Xriyi = Xr2,r'2 = Xrr' (15) 
Mrl = A*r2 = (16) 
Xrl,r2^Xv (17) 

Here, Xv is real and positive, and all other inter-layer x's 



7 




(b) 

FIG. 3. Unit cells used for SCM calculations on the cubic 
lattice (a), single bilayer (b), and single-layer square lattice 
(c). In the cubic case the primitive Bravais lattice vectors 
are chosen parallel to the edges of the rectangular prismatic 
unit cell. The analogous statement is true for the bilayer 
and single-layer cases, with primitive Bravais lattice vectors 
parallel to the l^^y edges of the unit cell. 



are assumed to vanish. We let n label the one-particle 
eigenstatcs of a single layer, with energies . The full 
one-particle spectrum is then given by 



,2d 



(18) 



where cr = ±1. We assume that the energy spectrum 
and filling are such that only a = —1 states are oc- 
cupied by fermions, in which case the two-dimensional 
spectrum (shifted in energy by —Xv) is filled by 
NNs/k = 2NN'^'^/k fermions. This corresponds to a 
single-layer problem with twice as many fermions, or fill- 
ing parameter k' — k/2. The saddle point energy is then 



Em FT = NN. 

2N 



2 

2d 

J 

{rr') 



-Xv 



(19) 



Ef{k') 



Here, 



= 2771, and E'j^{k') is the ground state energy 
of the fermionic part of the mean-field Hamiltonian [last 
two terms of Eq. (8)] , for a single-layer square lattice with 
filling parameter k'. The first two terms of Eq. (19) are 
minimized with respect to Xv to find Xv = J/k. The 
last three terms combine to Elfpj,{k',J''), the saddle 
point energy of a single-layer square lattice with filling 



K 


Cubic lattice 


Single bilayer 


k/2 square lattice 


5 


-1 ^ ^x,y.z 5 


10 

30 


1 < lx,y < 5 


10 
60 


1 < lx,y < 6 


30 
60 


6 


1 ^ ^x.y,z ^ 6 


4 

30 


1 < lx,y < 6 


4 

60 






7 


1 ^ Ix.y ,z ^ 7 


4 
21 


1 < L,y < 7 


10 

35 


l<lx<7 

1 < ly < 10 


20 

42 


8 


1 *^ ^x,y ^ 8 


4 
24 


1 < lx,y < 8 


4 

40 






9 


1 < l^,y < 9 
l<h<4: 


4 

36 


1 < < 9 

1 < ly < 11 


10 

36 


1 < < 10 

l<ly <9 


10 

36 


10 


1 < Ix.y < 10 

1<IZ<4: 


4 

30 


1 < L,y < 10 


5 

60 







TABLE II. This table contains information about our SCM 
numerical study on the cubic lattice (1st column), as well as 
the related problems of a single bilayer (2nd column), and 
single layer square lattice with k' = k/2 (3rd column). On 
the left-hand side of each entry of the table, the range of unit 
cell dimensions is shown as an inequality. For every choice of 
lx,y,z within the given range, the number of times we ran the 
SCM algorithm with distinct random initial configurations of 
Xrr' is shown on the right-hand side of the entry (top). Also 
on the right-hand side is the minimum linear system size L 
(bottom, italics). 



parameter k' and J"' ~ J /2. Noting that 



(20) 



we obtain the following relation between bilayer and 
single-layer saddle point energies: 



Emft 
N,N 



3_ 

2fc2 



i i;|,V(fc/2) 

4 N1<^N 



(21) 



This relation allows us to study via SCM the single- 
layer square lattice with filling parameter k' = k/2 as 
a further check on the cubic lattice results. For inte- 
ger fc', this was already done in Ref. 7. We carried out 
SCM calculations for the half-odd integer filling parame- 



ters k' 



(see Table II and Fig. 3). For all values 



of fc, these calculations find the same ground states as 
found by the single-bilayer SCM calculations. 

As a further check on our results, we also computed 
the energies of some simple competing states. Table III 
compares the energies of these states to the ground state 
saddle point energies found by SCM. 



IV. DISCUSSION 

The large- iV results presented here find a rich variety 
of candidate non-magnetic ground states for Mott insula- 
tors of ultra-cold fermionic AEA. It would be fascinating 
to realize any of these states experimentally. In order to 
achieve this, there still need to be substantial advances in 
preparation of low-entropy magnetic states of ultra-cold 
atoms, and our results add to the increasing motivation 
to pursue such advances specifically in AEA systems. In 
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K 
O 


u 


7 


B 
o 




1 n 


Bilayer (<3> — 2nn/k) 


-0.0444916 


-0.0344012 


-0.0273888 


-0.0223613 


-0.0186271 


-0.01577 


fc-site cluster 


-0.04 


-0.032407 


-0.026239 


-0.0234375 


-0.0178326 


-0.014 


Uniform real x 


-0.0394159 


-0.0312776 


-0.0254048 


-0.0210391 


-0.0177088 


-0.0151133 




-0.0430802 


-0.0330693 


-0.0261299 


-0.0212772 


-0.0177579 


-0.0151134 


SCM ground state 


-0.0445021 


-0.0347222 


-0.0273888 


-0.0234375 


-0.0188265 


-0.01577 



TABLE III. Comparison of energies of a variety of simple saddle points (top four rows), with the energy of the ground state found 
by SCM numerics (bottom row). All energies are in units of NJ'Ns- Each row represents a class of saddle points, described 
below. For classes including multiple different saddle points, the energy shown is the lowest in the class. We considered the 
following classes of saddle points: Bilayer ($ — 2nn/k). We considered a generalization of the CSL bilayer saddle point 
described in the main text, where the flux through each plaquette is "I> = 2im/k, where n = 0,...,fc — 1. k-site cluster. The 
energy of a cluster with k sites is proportional to the number of bonds in the cluster,''^ ' so the lowest-energy such state can be 
found by finding a fc-site cluster containing the greatest number of bonds. Uniform real x- This is the state where Xrr' is real 
and spatially constant. {2Tinx,y,z/k)-flux. These states have 2'Knx/k flux through every plaquette normal to the i-direction, 
and similarly for y and z, where < nx,y,z <k — l. Since most of these states break lattice rotation symmetry, the magnitude 
IXrr'l is allowed to vary depending on bond orientation, but is fixed to be translation invariant.^" 



addition, if future experiments ean enter a regime where 
any of the states diseussed here can be realized, it will 
be of crucial importance to devise probes of their char- 
acteristic properties. 

We would like to close by highlighting the CSL bilayer 
state, which has some striking properties that would be 
fascinating to realize experimentally, and which we now 
briefly discuss. At the large- mean-field level the cu- 
bic lattice breaks into disconnected bilayers, and one can 
understand the properties beyond mean-field theory by 
first focusing on a single bilayer. The effect of fluctua- 
tions is to couple the fermions to a dynamical U(l) gauge 
field. The mean-field fermions are in a gapped integer 
quantum Hall state, so integrating them out generates 
a Chcrn-Simons term for the U(l) gauge field. Because 
the mean-field fermions in a single bilayer and in the 
single-layer square lattice CSL'' '^ ' have in both cases a 
single chiral edge mode per spin species, the coefficient of 
the Chern-Simons term and associated topological prop- 
erties are the same. The spinous are Abelian anyons 
with statistics angle = tt ± ir/N , and there is a chiral 
edge mode with gapless excitations carrying SU{N) spin, 
which is described by a chiral SU{N)i Wcss-Zumino- 
Witten model.''' ' 



If adjacent bilayers are coupled weakly, bulk properties 
are unaffected due to the energy gap. One simply has a 
many-layer CSL state, with anyonic spinous confined to 
the the individual layers. Due to the gapless edge modes 
of single bilayers, the physics on the two-dimensional sur- 
face is likely more interesting. This depends crucially 
on whether adjacent bilayers have the same or opposite 
magnetic fiux, as the direction of the fiux controls the di- 
rection of the chiral edge modes. If the fluxes are aligned 
oppositely in neighboring bilayers, then edge modes on 
neighboring bilayers are counterpropagating and an en- 
ergy gap is possible on the two-dimensional surface. On 
the other hand, if all fluxes are parallel, then all the chi- 
ral edge modes propagate in the same direction, and the 
two-dimensional surface is expected to remain gapless. 
The resulting surface state is a kind of two-dimensional 
chiral "spin metal," which could be interesting to study 
in future work. 
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